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Entanglement renormalization techniques are applied to numerically investigate the ground state 
of the spin- 1 Heisenberg model on a kagome lattice. Lattices of = {36, 144, oo} sites with periodic 
boundary conditions are considered. For the infinite lattice, the best approximation to the ground 
state is found to be a valence bond crystal (VBC) with a 36-site unit cell, compatible with a previous 
proposal. Its energy per site, E = —0.43221, is an exact upper bound and is lower than the energy 
of any previous (gapped or algebraic) spin liquid candidate for the ground state. 
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Low dimensional spin-^ quantum systems have long 
been the focus of intense research efforts, largely fueled 
by the search for exotic states of matter. An important 
example of a geometrically frustrated quantum antiferro- 
magnet |lj is the spin- ^ kagome-lattice Heisenberg model 
(KLHM). Despite a long history of study, the nature of 
its ground state remains an open question. Leading pro- 
posals include valence bond crystal (VBC) [2-7 and spin 
liquid (SL) [8:16 ground states. Interest has been fur- 
ther stimulated by recent experimental work on Herbert- 
smithite ZnCu3(OH)6Cl2, a possible physical realization 
of the model [17 . 

Progress in our understanding of the KLHM has been 
hindered, as with many other models of frustrated anti- 
ferromagnets, by the inapplicability of quantum Monte 
Carlo methods due to the negative sign problem. Never- 
theless, systems with up to 36 sites have been addressed 
with exact diagonalization [9", T8 , whereas the density 
matrix renormalization group (DMRG) has been used 
to explore lattices of order N ^ 100 sites [16 . Unfortu- 
nately, it is very difficult to infer the nature of the ground 
state of an infinite system from these results. The rea- 
son is that these lattices are still relatively small given 
the 36-site unit cell of the leading VBC proposal, or the 
algebraic decay of correlations in some SL proposals. In 
larger systems, support for a SL ground state has also 
been obtained with a SL ansatz p!3HT5] , whereas evidence 
for a VBC has been obtained for an infinite lattice with 
a series expansion around an arbitrary dimer covering 
iniEs, but both approaches are clearly biased. 

In this paper we report new, independent numerical 
evidence in favor of a VBC ground state for the KLHM 
model. This is done with entanglement renormalization 
[T9H25], a real space RG approach that, through the 
proper removal of short-range entanglement, is capable 
of providing an approximation to ground states of large 
2D lattices [24| by means of the multi-scale entanglement 
renormalization ansatz (MERA) [20] . After describing a 
scheme for the Kagome lattice with periodic boundary 
conditions, we address lattices of N = {36, 144, oo} sites. 
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FIG. 1: (Color on-line) Coarse-graining transformation that 
maps (i) a kagome lattice Co of N sites into (vi) a coarser 
lattice Ci of A^/36 sites, (ii) The lattice is first partitioned 
into blocks of 36 sites, (iii) Disentanglers u are applied across 
the corners of three blocks, followed by (iv) disentanglers v 
applied across the sides of two neighboring blocks, (v) Isome- 
trics w map blocks to an effective site of the coarse-grained 
lattice, (vii) Tensors u, v and w have a varying number of 
incoming and outgoing indices according to Eq. [l] 



Our simulations converge to a VBC state compatible with 
that first proposed by Marston and Zeng [2^ and revis- 
ited by Nikolic and Senthil [4 , and by Singh and Huse 
[6l [7] . For an infinite lattice we obtain an energy per site 
E = —0.43221. This energy corresponds to an explicit 
(MERA) wave-function and therefore provides us with 
a strict upper bound for the true ground state energy. 
Importantly, its value is lower than the energy of any ex- 
isting SL ansatz on a sufficiently large lattice, which we 
interpret as strong evidence for a VBC ground state in 
the thermodynamic limit. Our results are also the first 
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demonstration of the utility of entanglement renormal- 
ization to study 2D lattice models that are beyond the 
reach of quantum Monte Carlo techniques. 

The present approach is based on the coarse-graining 
transformation of Fig. [l] which is applied to a kagome 
lattice Cq made of N sites. It maps blocks of 36 sites 
of Cq onto single sites of a coarser lattice Ci made of 
sites. A Hamiltonian Hq defined on lattice Cq be- 
comes an effective Hamiltonian Hi on lattice Ci. Anal- 
ogously, the ground state |^o) of is transformed into 
the ground state of i^i. The transformation decom- 
poses into three steps. Firstly disentanglers unitary 
tensors that act on 9 sites, are applied across the cor- 
ners of three neighboring blocks. Then disentanglers v 
are applied across the sides of two neighboring blocks; 
these tensors reduce ten sites (each described by a vector 
space C2 of dimension 2) into two effective sites (each 
described by a vector space of dimension x)- Finally 
isometrics w map the remaining sites of each block into 
a single effective site of Ci. Thus the tensors li, v and 
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transform the ground state |^o) of lattice Cq into the 
ground state of lattice Ci through the sequence 

1*0) 4 1*;,) 4 1*;,') ^ I*,). (2) 

The disentanglers u and v aim at removing short-range 
entanglement across the boundaries of the blocks; there- 
fore states I^q) and |^q) possess decreasing amounts of 
short range entanglement. If state |^o) only has short- 
range entanglement to begin with, then it is conceivable 
that the state has no entanglement left at all. For 
a finite lattice {N = 144) we consider a state |^o) that, 
after the coarse-graining transformation, gives rise to an 
entangled state on = 4 sites. For an infi- 

nite lattice we make instead an important assumption, 
namely that is a product (non-entangled) state (this 
is equivalent to setting x = 1 in Eq. [T]). How short- 
ranged must the entanglement in |^o) be for this as- 
sumption to be valid? By reversing the transformation 
on a product state it can be seen that each site in 

1^0) is still entangled with at least 84 neighboring sites. 

The disentanglers and isometrics (u^v^w) were initial- 
ized randomly and then optimized so as to minimize the 
expectation value of the KLHM Hamiltonian, 

Hq = J ^ Si • Sj, (3) 

by following the algorithms of Ref. [25l, with cost 
0{2^'^X^X^) [28 . Specifically, for lattices with N = 36 
and 144 sites, the resulting (one-site and four-site) Hamil- 
tonian Hi is diagonalized exactly. Instead, for = 00, 




FIG. 2: (Color on-line) The 36-site unit cell for the honey- 
comb VBC, strong bonds are drawn with thick lines. Three 
different types of strong bonds can be identified; the six bonds 
belonging to the pinwheels (red) , six bonds belonging to each 
'perfect hexagon' (blue) and the parallel bonds between per- 
fect hexagons (green) . Dotted arrows indicate the axis where 
spin- spin correlators have been computed. Bond-bond cor- 
relators have been computed between the reference bond (1) 
and the other numbered bonds. 



we use the finite correlation range algorithm (Sect. V.D 
of ref. [25 ). All computations led to highly dimerized 
wave-functions of the VBC type. In order to explain the 
results, consider the exact 'honeycomb' VBC state, de- 
noted |h-VBC), whose 36-site unit cell is shown in Fig. [2j 
Each unit cell contains two 'perfect hexagons' (resonat- 
ing bonds around a hexagon) and a 'pinwheel'. Three 
different types of strong bonds can be identified: those 
of the pinwheels (red), parallel bonds (green) and per- 
fect hexagons (blue). The pinwheel and parallel bonds 
are singlets (energy per bond = —0.75) while the perfect 
hexagons are in the ground state of a periodic Heisenberg 
chain of 6 sites (energy per bond = —0.4671). The rest 
of links have zero energy. We call a 'honeycomb' VBC 
a state that has strong bonds according to the above 
pattern, even though the rest of bonds (weak bonds) 
need not have zero energy. The 'honeycomb' VBC was 
originally proposed by Marston and Zeng [2 (see also 
HEl C]). Our simulations with N = 144 and A' = 00 
produce a VBC of this type as the best MERA approxi- 
mation to the ground state. 

The energies obtained for an infinite lattice are shown 
in Table [l| For each value of x (see Eq. [l]), the MERA is 
an explicit wave-function and therefore provides an up- 
per bound to the exact ground state energy. Energies 
computed for the A^ = 144 lattice matched those of the 
infinite lattice to within 0.02% and have been omitted. 
Our N = 00 energy is also compatible with the energy 
E = -0.433 ± 0.001 for a VBC obtained with series ex- 
pansion in Ref. |6|, and is lower than that obtained in 
Ref. with DMRG {E = -0.43160 for A" = 108) 
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FIG. 3: (top) Bond energies for the 36-site unit cell of infi- 
nite MERA wave- functions, for two different values of x, as 
compared to those of an exact honeycomb VBC, |h-VBC), 
and those of a spin liquid, which by definition has all equal 
strength bonds. The MERA wave-functions clearly match 
the proposed honeycomb VBC; we identify (i) the six strong 
'pinwheel' bonds (red bonds), the six 'parallel' bonds (green 
bonds) and (iii) the 12 'perfect hexagon' bonds (blue bonds). 
The (iv) remaining 48 bonds are the weak bonds of the 
unit cell, (bottom) Bond energies for the 36-site lattice. 
Here a randomly initialized MERA converges to a dimerized 
state that does not match the honeycomb VBC pattern, but 
gives lower overall energy than a honeycomb VBC initialized 
MERA of the same x- 



and in Ref. [14 with fermionic mean- field theory and 
Gutzwiller projection {E = -0.42863 for N = 432). We 
further notice that, where finite size effects are still rele- 
vant, such as in the N = 108 case, they tend to decrease 
the ground state energy. 

Fig. [3]shows the distribution of bond energies obtained 



TABLE I: Ground state energies as a function of x- 
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-0.42145 


-0.42164 


-0.42143 
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-0.42952 
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FIG. 4: Spin-Spin correlators along arrows 'A' and 'B' of Fig. 
|2]for infinite lattice MERA of x = 4 and x = 16. Although 
along both lattice directions considered the correlators decay 
exponentially, the decay along arrow 'A' (a line joining two 
perfect hexagons) is seen to be slower than along arrow 'B' 
(a line joining perfect hexagon to pinwheel). The plateaus 
marked (i), (ii) and (iii) show the correlation is the same with 
both spins of a strong bond. 



for the N = oo lattice. With x = 4, one observes an en- 
ergy increase per site over |h-VBC) of 0.08 in the par- 
allel (green) bonds and also in some of the hexagon (blue) 
bonds, with the weak bonds having lower energy in re- 
turn. As X is increased, the energy of the 'strong' bonds 
becomes slightly larger and that of the 'weak' bonds con- 
tinues to decrease. However, the dimerization clearly sur- 
vives: the bond energies are not seen to converge to a 
uniform distribution as required for a SL. 

Fig. |4] shows spin-spin correlators evaluated along two 
different lattice axis A and B (cf. Fig. |2| for A/" = oo. 
These correlators decay exponentially with well defined 
'plateaus', where the correlation is the same with both 
spins of a strong bond. Correlations along the line join- 
ing a perfect hexagon and a pinwheel are seen to decay 
faster than along the line joining two perfect hexagons, 
consistent with the observation from Fig. |3]that the pin- 
wheel bonds remain almost exact singlets even for high 
values of x- Table [ll] shows bond-bond connected and 
disconnected correlators, Ci^a and Di ^^ 

C,. . {{S.S)^{S.S)J (4) 

Di,a = Ci,„-((5-5)^)((5-5)J, (5) 

between a reference bond '1' and a surrounding bond 
a = 1, • • • ,14 (cf. Figj2|. While disconnected correlators 
decay exponentially with distance, some connected cor- 
relators remain significant at arbitrary distances, demon- 
strating the long-range order of the VBC state. 

Let us discuss the results for a lattice with = 36 
sites. When initialized with random tensors, the MERA 
produced VBC type configurations which typically did 
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TABLE II: Bond-Bond correlators for x = 16 (A^ = oo). 



Bond 








1 


0.38639 


0.22815 


-0.39780 


2 


0.29490 


0.03496 


-0.65342 


3 


0.00422 


-0.00115 


-0.01349 


4 


0.09503 


0.04644 


-0.12211 


5 


0.22488 


0.06609 


-0.39918 


6 


-0.00609 


-0.01426 


-0.02054 


7 


-0.10372 


-0.03376 


-0.34561 


8 


0.26083 


0.01252 


-0.62421 


9 


0.00336 


0.00368 


0.00082 


10 


-0.00001 


-0.00017 


-0.00042 


11 


0.29113 


0.00196 


-0.72693 


12 


0.00090 


-0.00036 


-0.00321 


13 


0.16632 


0.00026 


-0.41744 


14 


0.15674 


-0.00037 


-0.39496 



not match the honeycomb VBC, although simulations 
initialized in the state h-VBC) retained the honeycomb 
VBC pattern (see Fig. 3|. Here the randomly initialized 
VBC produced a lower energy (0.5% above the exact di- 
agonalization result E = —0.438377 of [9 ) than the hon- 
eycomb VBC type solution for an equivalent value of x 
(cf. Table |l]) . These results strongly suggest that finite 
size effects in the N = 36 site lattice lead to a significant 
departure from the physics of the infinite system. 

To summarize, we have used entanglement renormal- 
ization techniques to obtain new numerical evidence in- 
dicating that the ground state of the KLHM is of the 
honeycomb VBC type. In order to assess the robustness 
of this result, we briefly discuss some of the limitations 
of the present approach. 

Firstly, the coarse-graining transformation of Fig. [T] 
which maps 36 sites into one site, was designed to en- 
sure compatibility with the 36-site unit cell of honey- 
comb VBC type solutions. While our approach did not 
preclude solutions with a smaller, compatible unit cell 
(such as a 12-site unit cell or a fully translation invariant 
solution), we cannot rule out the possibility that a state 
with an incompatible unit cell might have a lower energy. 

Secondly, the infinite lattice was investigated by re- 
stricting the range of entanglement in the ansatz to 
blocks of 84 spins, imposed through an unentangled state 
l^i) in Eq. [2] This restriction was only implemented af- 
ter preliminary simulations with x = 12 had produced 
identical energies irrespectively of whether was al- 
lowed to be entangled. However, it could still be that 
entanglement in would make a big difference for 

larger values of x- We find this scenario quite unlikely, 
but could not test it due to computational limitations. 

Finally, the MERA is an essentially unbiased method 
provided that the candidates to be the ground state of 



the system have all a relatively small amount of entangle- 
ment. But when deciding between a VBC (which mostly 
has short-range entanglement) and e.g. the algebraic 
SL of Refs. [141 US] (significantly more entangled at all 
length scales), it might well be that the MERA is biased 
toward the low entanglement solution. Therefore our re- 
sults do not conclusively exclude a SL ground state. We 
emphasize, however, that the ground state energies ob- 
tained with the MERA are lower than the SL energies of 
Refs. [HHTS]. 

We thank P. Lee, R. Singh, X.-G. Wen and S. White 
for valuable conversations. Support from the Australian 
Research Council (APA, FF0668731, DP0878830) is ac- 
knowledged. 
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Appendix A: The Triangular Lattice 

In this appendix we describe a coarse-graining trans- 
formation for the triangular lattice that, when used in 
conjunction with the scheme of Fig. [l] allows one to an- 
alyze arbitrarily large or infinite kagome lattices without 
restricting in Eq. |2]to be unentangled. This scheme 
may also be used directly to simulate models defined on 
a triangular lattice. 



(i) Lattice, 




FIG. 5: (Color on-line) A coarse-graining scheme which trans- 
forms (i) a triangular lattice £r of N sites into (vi) a coarser 
triangular lattice £t+i of sites, (ii) Firstly the lattice 

is partitioned into blocks of 16 sites, (iii) Disentanglers u 
are applied across the corners of three blocks, followed by 
(iv) disentanglers v applied across the sides of two neighbor- 
ing blocks, (v) Isometrics w are applied to map blocks to 
an effective site of the coarse-grained lattice, (vii) Tensors 
V and w have a varying number of incoming and outgoing 
indices according to Eq. [6] 



Fig. [5] shows a coarse-graining transformation that 
maps a triangular lattice Cr of N sites into and a coarser 
triangular lattice jCr-\-i of N/16 sites. This transforma- 
tion is accomplished in three steps; firstly disentanglers 

unitary tensors which act on three sites, are applied 
across the corner of three blocks. Next disentanglers v 
are applied across the sides of neighboring blocks; these 
are isometric tensors that truncate five lattice sites (each 
one described by a vector space of dimension x) iiito 
two effective sites (each one described by a vector space 
C-^ of dimension x). Finally, isometrics w map the re- 
maining sites of each block into a single effective site of 
Cr-\-i- Collectively, the (isometric) tensors v and 



^ ^x 
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(6) 



transform the state 1^^-) G C^*^ of the lattice Cr into a 
state l^r+i) e C^®^^^^ of the effective lattice £^+1. 



A MERA for the state |^o) of kagome lattice system 
can be built by composing the transformation of Fig. [l] 
with several transformations as in Fig. |5j Following this 
recipe, c.f. Fig. |6j a MERA of r layers allows for investi- 
gations of finite kagome lattices of A/" = 36 x 16'^""'^ sites. 
This scheme can also be used with the scale-invariant 
algorithm (Sect V.C of [^) to investigate scale- invariant 
systems, such as critical systems, on infinite kagome 
lattices. The triangular coarse-graining transformation 
of Fig. [5] may also be used directly to investigate models 
defined on triangular lattices. A MERA with r layers, 
each defined by the transformation Fig. [5) represents a 
wave-function on a triangular lattice of = 16^^ sites. 
Indeed, when applied to a spin-^ Heisenberg model on 
an TV = 00 triangular lattice, this scheme reproduces 
the three-sublattice 120° order found in other studies 
and gives a ground energy for the infinite system at 
E = —0.545 per site, comparable to a series expansion 
estimate oi E = -0.5502 per site 127]. 
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FIG. 6: (i) A kagome lattice MERA can be built using the 
transformation of Fig. [l] as the bottom MERA layer, equiv- 
alent ly as the initial coarse- graining transformation, and the 
transformation of Fig. [s] for subsequent MERA layers, (ii) 
The transformation of Fig. |5] can also be used directly to 
investigate models defined on triangular lattice geometries. 



Appendix B: Defective Crystals 

This appendix details a modification of the MERA op- 
timization algorithm [25] that was found necessary to 
ensure proper convergence of some simulations of the 
KLHM. 
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N = 144, defective valence bond cystal 




N = 144, honeycomb valence bond cystal 




FIG. 7: (Color on-line) Converged x = 8 MERA solutions for 
the N = 144 site KLHM. The thickness of a line between two 
lattice sites is proportional to the absolute bond energy, \E\. 
Blue solid lines are used for negative energy bonds, E < 0, 
and red dashed lines for positive energy bonds, ^ > 0. (Top) 
An example of a 'defective' VBC. Local structures that match 
the 'pinwheels' and 'perfect hexagons' of the honeycomb VBC 
are observed, but long range order is absent. (Bottom) By re- 
stricting the rate of lattice 'cooling', as per Eq|7| one more 
reliably obtains the honeycomb VBC from a randomly initial- 
ized MERA. On average a defective VBC had 0.05% greater 
energy than the honeycomb VBC of the same value of x- 



All simulations of the KLHM converged to highly 
dimerized VBC type states. For lattices of = {144, oo} 
sites the lowest energy state that was found, or best ap- 
proximation to the ground-state, had a pattern of strong 
bonds matching the 'honeycomb' VBC state. However, 
not all randomly initialized simulations converged to this 
state; in some instances the MERA wave- function con- 
verged to a 'defective' VBC state, an example of which 
is shown Fig. [7[ The observed defective VBC's ranged 
from states which matched the honeycomb VBC with 
only a few misplaced strong bonds to states with an al- 
most completely disordered placement of strong bonds. 
A defective VBC was observed to have on average 0.05% 
higher energy per site than the corresponding honeycomb 
VBC type MERA of the same x- The exact difference 



in energy depended on the amount of defects; typically, 
states that deviated more from the honeycomb VBC had 
higher energy than states which deviated less. In the 
instances that the method did not converge to the hon- 
eycomb VBC one may conclude that the optimization of 
the MERA became trapped in a local minimum. 

Becoming trapped in a local minimum is not surpris- 
ing given that the optimization of the MERA is based 
upon updates of individual tensors. During the initial 
iterations of the optimization, locally stable structures, 
such as singlets between neighboring spins, form. But 
they do not necessarily form in places compatible with 
a global minimization of the energy. Once the method 
has converged to a particular VBC state, the transfor- 
mation required to bring the state to a different VBC, 
with a new pattern of strong bonds, requires simultane- 
ous shifts of many strong bonds. This is unlikely to occur 
with an optimization based on energy- lowering updates 
of individual tensors. 

There are many ways to decrease the risk of becom- 
ing trapped in a local minimum due to the formation of 
local singlets. In the present work this is achieved by re- 
stricting how much a tensor is allowed to change in one 
single update; this may be thought of as decreasing the 
rate at which the MERA is 'cooled-down' from an initial 
(high-energy) state to the (low-energy) approximation to 
the ground-state. The goal of this restriction is to pre- 
vent locally stable structures from forming too early in 
the optimization. In order to discuss how this restric- 
tion may be implemented, we first refresh some details of 
the variational MERA optimization, see Sect. IV of [25] . 
In the standard algorithm, a tensor w in the MERA is 
replaced with a better (that is leading to lower energy) 
tensor w' as follows. First one computes the linearized 
environment of denoted T^. Then, given the singular 
value decomposition of the environment, T^^ = USV^^ 
one chooses the updated tensor to be = VU^ . Here 
we wish to restrict the amount e the tensor can change, 
lliL' — < e. This is achieved by using a modified en- 
vironment defined as 

= + (7) 

for some A > 0. As before, given the singular value de- 
composition of the modified environment = USV^, 
the updated tensor is chosen to be = VU^ . The 
parameter A acts as a soft constraint to ensure the tensor 
change e is kept small with each update, a larger value of 
A giving a smaller change e. In the limit of A ^ oc, the 
tensor w would remain constant, e = 0, in the update. 
Simulations of the KLHM with a suitable value of A, say 
A > 10 for the initial optimization, were much more likely 
to converge to a defect free honeycomb VBC than were 
unrestrained (A = 0) optimizations. 



